{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IYul_lq4tn3D"
   },
   "source": [
    "##### Copyright 2020 The OpenFermion Developers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "id": "sP7GfWPstpBt"
   },
   "outputs": [],
   "source": [
    "#@title Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "#\n",
    "# https://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "a43053f38ce0"
   },
   "source": [
    "# Introduction to the bosonic operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8R3_GRhvt88K"
   },
   "source": [
    "<table class=\"tfo-notebook-buttons\" align=\"left\">\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://quantumai.google/openfermion/tutorials/bosonic_operators\"><img src=\"https://quantumai.google/site-assets/images/buttons/quantumai_logo_1x.png\" />View on QuantumAI</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://colab.research.google.com/github/quantumlib/OpenFermion/blob/master/docs/tutorials/bosonic_operators.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/colab_logo_1x.png\" />Run in Google Colab</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://github.com/quantumlib/OpenFermion/blob/master/docs/tutorials/bosonic_operators.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/github_logo_1x.png\" />View source on GitHub</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a href=\"https://storage.googleapis.com/tensorflow_docs/OpenFermion/docs/tutorials/bosonic_operators.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/download_icon_1x.png\" />Download notebook</a>\n",
    "  </td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "n5kB9jHTuI6B"
   },
   "source": [
    "## Setup\n",
    "\n",
    "Install the OpenFermion package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "5cbe6b680387"
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import openfermion\n",
    "except ImportError:\n",
    "    !pip install git+https://github.com/quantumlib/OpenFermion.git@master#egg=openfermion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "952ae604e015"
   },
   "source": [
    "## The BosonOperator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "e41da94a83e0"
   },
   "source": [
    "Bosonic systems, like Fermionic systems, are expressed using the bosonic creation and annihilation operators $b^\\dagger_k$ and $b_k$ respectively. Unlike fermions, however, which satisfy the Pauli exclusion principle and thus are distinguished by the canonical fermionic anticommutation relations, the bosonic ladder operators instead satisfy a set of commutation relations:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& [b_i^\\dagger, b_j^\\dagger] = 0, ~~~ [b_i, b_j] = 0, ~~~ [b_i, b^\\dagger_j] = \\delta_{ij}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Any weighted sums of products of these operators are represented with the `BosonOperator` data structure in OpenFermion. Similarly to when we introduced the `FermionOperator`, the following are examples of valid `BosonOperator`s:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& a_1 \\nonumber \\\\\n",
    "& 1.7 b^\\dagger_3 \\nonumber \\\\\n",
    "&-1.7 \\, b^\\dagger_3 b_1 \\nonumber \\\\\n",
    "&(1 + 2i) \\, b^\\dagger_3 b^\\dagger_4 b_1 b_9 \\nonumber \\\\\n",
    "&(1 + 2i) \\, b^\\dagger_3 b^\\dagger_4 b_1 b_9 - 1.7 \\, b^\\dagger_3 b_1 \\nonumber\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The `BosonOperator` class is contained in `ops/_boson_operators.py`. The `BosonOperator` is derived from the `SymbolicOperator`, the same class that derives the `FermionOperator`. As such, the details of the class implementation are identical - as in the fermion case, the class is implemented as hash table (python dictionary). The keys of the dictionary encode the strings of ladder operators and values of the dictionary store the coefficients - the strings are subsequently encoded as a tuple of 2-tuples which we refer to as the \"terms tuple\".\n",
    "\n",
    "\n",
    "Each ladder operator is represented by a 2-tuple. The first element of the 2-tuple is an int indicating the quantum mode on which the ladder operator acts. The second element of the 2-tuple is Boole: 1 represents raising and 0 represents lowering. For instance, $b^\\dagger_8$ is represented in a 2-tuple as $(8, 1)$. Note that indices start at 0 and the identity operator is an empty list. \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "I & \\mapsto () \\nonumber \\\\\n",
    "b_1 & \\mapsto ((1, 0),) \\nonumber \\\\\n",
    "b^\\dagger_3 & \\mapsto ((3, 1),) \\nonumber \\\\\n",
    "b^\\dagger_3 b_1 & \\mapsto ((3, 1), (1, 0)) \\nonumber \\\\\n",
    "b^\\dagger_3 b^\\dagger_4 b_1  b_9 & \\mapsto ((3, 1), (4, 1),  (1, 0), (9, 0)) \\nonumber\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Alternatively, the `BosonOperator` supports the string-based syntax introduced in the `FermionOperator`; in this case, the terms are separated by spaces, with the integer corresponding to the quantum mode the operator acts on, and `'^'` indicating the Hermitian conjugate:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "I & \\mapsto \\textrm{\"\"} \\nonumber \\\\\n",
    "b_1 & \\mapsto \\textrm{\"1\"} \\nonumber \\\\\n",
    "b^\\dagger_3 & \\mapsto \\textrm{\"3^\"} \\nonumber \\\\\n",
    "b^\\dagger_3 b_1 & \\mapsto \\textrm{\"3^}\\;\\textrm{1\"} \\nonumber \\\\\n",
    "b^\\dagger_3 b^\\dagger_4  b_1 b_9  & \\mapsto \\textrm{\"3^}\\;\\textrm{4^}\\;\\textrm{1}\\;\\textrm{9\"} \\nonumber\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "   Note that, unlike the `FermionOperator`, the bosonic creation operators of different indices commute. As a result, the `BosonOperator` automatically sorts groups of annihilation and creation operators in ascending order of the modes they act on.\n",
    "</div>\n",
    "\n",
    "\n",
    "Let's initialize our first term! We do it two different ways below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "39d1ec24c9d5"
   },
   "outputs": [],
   "source": [
    "from openfermion.ops import BosonOperator\n",
    "\n",
    "my_term = BosonOperator(((3, 1), (5, 0), (4, 1), (1, 0)))\n",
    "print(my_term)\n",
    "\n",
    "my_term = BosonOperator('3^ 5 4^ 1')\n",
    "print(my_term)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0a9e12ccd92d"
   },
   "source": [
    "Note the printed order differs from the code, since bosonic operators of different indices commute past each other.\n",
    "\n",
    "The preferred way to specify the coefficient in openfermion is to provide an optional coefficient argument. If not provided, the coefficient defaults to 1. In the code below, the first method is preferred. The multiplication in the second method actually creates a copy of the term, which introduces some additional cost. All inplace operands (such as +=) modify classes whereas binary operands such as + create copies.\n",
    "\n",
    "The additive and multiplicative identities can also be created:\n",
    "\n",
    "* `BosonOperator(())` and `BosonOperator('')` initialises the identity (`BosonOperator.identity()`).\n",
    "* `BosonOperator()` and `BosonOperator()` initialises the zero operator (`BosonOperator.zero()`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "dac42b48422e"
   },
   "outputs": [],
   "source": [
    "good_way_to_initialize = BosonOperator('3^ 1', -1.7)\n",
    "print(good_way_to_initialize)\n",
    "\n",
    "bad_way_to_initialize = -1.7 * BosonOperator('3^ 1')\n",
    "print(bad_way_to_initialize)\n",
    "\n",
    "identity = BosonOperator('')\n",
    "print(identity == BosonOperator.identity())\n",
    "print(identity)\n",
    "\n",
    "zero_operator = BosonOperator()\n",
    "print(zero_operator == BosonOperator.zero())\n",
    "print(zero_operator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "97987980c581"
   },
   "source": [
    "Note that `BosonOperator` has only one attribute: .terms. This attribute is the dictionary which stores the term tuples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "d0fdd98b0012"
   },
   "outputs": [],
   "source": [
    "my_operator = BosonOperator('4^ 1^ 3 9', 1. + 2.j)\n",
    "print(my_operator)\n",
    "print(my_operator.terms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dd0f0b072be5"
   },
   "source": [
    "### Methods and functions that act on the BosonOperator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "532f85ad65ab"
   },
   "source": [
    "There are various functions and methods that act on the `BosonOperator`; these include the ability to normal order, double check if the operator is Hermitian, and calculate the Hermitian conjugate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "52798c455b9b"
   },
   "outputs": [],
   "source": [
    "from openfermion.utils import hermitian_conjugated, is_hermitian\n",
    "from openfermion.transforms import normal_ordered"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "e5e898e7c471"
   },
   "source": [
    "`normal_ordered_boson` applies the bosonic commutation relations to write the operator using only normal-ordered terms; that is, that all creation operators are to the left of annihilation operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0339bdaf4842"
   },
   "outputs": [],
   "source": [
    "H = BosonOperator('0 0^', 1. + 2.j)\n",
    "H.is_normal_ordered()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "3d0f7bcb6a64"
   },
   "outputs": [],
   "source": [
    "normal_ordered(BosonOperator('0 0^', 1. + 2.j))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "69378cf1928a"
   },
   "source": [
    "We can also use a boson operator method to check if the operator conserves the particle number - that is, for each qumode, the number of annihilation operators equals the number of creation operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ef79ae4a3e1c"
   },
   "outputs": [],
   "source": [
    "H.is_boson_preserving()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0c6858ea0499"
   },
   "outputs": [],
   "source": [
    "H = BosonOperator('0 0^ 1^ ', 1. + 2.j)\n",
    "H.is_boson_preserving()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3dc0ec0da2e9"
   },
   "source": [
    "The Hermitian conjugated function returns the Hermitian conjugate of the operator, and its hermiticity can be checked using `is_hermitian`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "b9af0096bc62"
   },
   "outputs": [],
   "source": [
    "is_hermitian(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6554a0e92519"
   },
   "outputs": [],
   "source": [
    "hermitian_conjugated(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "946e2e2eafa5"
   },
   "outputs": [],
   "source": [
    "H = BosonOperator('0 1^', 1/2.)\n",
    "H += BosonOperator('1 0^', 1/2.)\n",
    "print(is_hermitian(H))\n",
    "print(hermitian_conjugated(H))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5421d9012891"
   },
   "source": [
    "## The QuadOperator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "d2e01ab8e1a2"
   },
   "source": [
    "Using the bosonic ladder operators, it is common to define the canonical position and momentum operators $\\hat{q}$ and $\\hat{p}$:\n",
    "\n",
    "$$ \\hat{q}_i = \\sqrt{\\frac{\\hbar}{2}}(\\hat{b}_i+\\hat{b}^\\dagger_i), ~~~ \\hat{p}_i = -i\\sqrt{\\frac{\\hbar}{2}}(\\hat{b}_i-\\hat{b}^\\dagger_i)$$\n",
    "\n",
    "These operators are Hermitian, and are referred to as the phase space quadrature operators. They satisfy the canonical commutation relation\n",
    "\n",
    "$$ [\\hat{q}_i, \\hat{p}_j] = \\delta_{ij}i\\hbar$$\n",
    "\n",
    "where the value of $\\hbar$ depends on convention, often taking values $\\hbar=0.5$, $1$, or $2$.\n",
    "\n",
    "In OpenFermion, the quadrature operators are represented by the `QuadOperator` class, and stored as a dictionary of tuples (as keys) and coefficients (as values). For example, the multi-mode quadrature operator $q_0 p_1 q_3$ is represented internally as `((0, 'q'), (1, 'p'), (3, 'q'))`. Alternatively, `QuadOperators` also support string input - using string input, the same operator is described by `'q0 p1 q3'`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "41892572d462"
   },
   "outputs": [],
   "source": [
    "from openfermion.ops import QuadOperator\n",
    "\n",
    "H = QuadOperator('q0 p1 q3')\n",
    "print(H)\n",
    "print(H.terms)\n",
    "\n",
    "H2 = QuadOperator('q3 p4', 3.17)\n",
    "H2 -= 77. * H\n",
    "print('')\n",
    "print(H2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bd9cf3d47089"
   },
   "source": [
    "Note that quadrature operators of different indices commute; as such, like the `BosonOperator`, by default we sort quadrature operators such that the operators acting on the lowest numbered mode appear to the left."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ab94957a0251"
   },
   "source": [
    "### Methods and functions that act on the QuadOperator\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "83ce3efcb3f8"
   },
   "source": [
    "Like the `BosonOperator`, there are various functions and methods that act on the `QuadOperator`; these include the ability to normal order, double check if the operator is Hermitian, and calculate the Hermitian conjugate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8a0b6b5a4b60"
   },
   "outputs": [],
   "source": [
    "from openfermion.utils import hermitian_conjugated, is_hermitian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "af10c754c9bc"
   },
   "source": [
    "`normal_ordered_quad` is an arbitrary convention chosen in OpenFermion that allows us to compare two quadrature operators that might be equivalent, but written in different forms. It is simply defined as a quadrature operator that has all of the position operators $\\hat{q}$ to the left of the momentum operators $\\hat{q}$. All quadrature operators can be placed in this 'normal form' by making use of the canonical commutation relation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bd5ead97318e"
   },
   "outputs": [],
   "source": [
    "H = QuadOperator('p0 q0', 1. + 2.j)\n",
    "H.is_normal_ordered()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cdca60a49d9b"
   },
   "outputs": [],
   "source": [
    "normal_ordered(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1e9731609ee8"
   },
   "source": [
    "By default, we assume the value $\\hbar=1$ in the canonical commutation relation, but this can be modified by passing the `hbar` keyword argument to the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "9e3d12b35141"
   },
   "outputs": [],
   "source": [
    "normal_ordered(H, hbar=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fe4a46aa5dad"
   },
   "source": [
    "We can also use a quad operator method to check if the operator is **Gaussian** - that is, all terms in the quad operator are of quadratic order or lower:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "0c4253372d9d"
   },
   "outputs": [],
   "source": [
    "H = QuadOperator('p0 q0', 1. + 2.j)\n",
    "H.is_gaussian()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "d2aa8581d1b1"
   },
   "outputs": [],
   "source": [
    "H = QuadOperator('p0 q0 q1', 1. + 2.j)\n",
    "H.is_gaussian()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6df5bae47a58"
   },
   "source": [
    "The Hermitian conjugated function returns the Hermitian conjugate of the operator, and its hermiticity can be checked using `is_hermitian`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "10029a18d3df"
   },
   "outputs": [],
   "source": [
    "H = QuadOperator('p0 q1 p1', 1-2j)\n",
    "hermitian_conjugated(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "954167974f66"
   },
   "outputs": [],
   "source": [
    "H = QuadOperator('p0 q0', 1/2.)\n",
    "H += QuadOperator('q0 p0', -1/2.)\n",
    "print(is_hermitian(H))\n",
    "print(hermitian_conjugated(H))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8fb0675ecb26"
   },
   "outputs": [],
   "source": [
    "H = QuadOperator('p0 q0', 1/2.)\n",
    "H += QuadOperator('q0 p0', 1/2.)\n",
    "print(is_hermitian(H))\n",
    "print(hermitian_conjugated(H))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "29e5c601ec1f"
   },
   "outputs": [],
   "source": [
    "hermitian_conjugated(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "50a84c0b5396"
   },
   "source": [
    "## Converting between quadrature operators and bosonic operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ae48bb4c0254"
   },
   "source": [
    "Converting between bosonic ladder operators and quadrature operators is simple - we just apply the definition of the $\\hat{q}$ and $\\hat{p}$ operators in terms of $\\hat{b}$ and $\\hat{b}^\\dagger$. Two functions are provided to do this automatically; `get_quad_operator` and `get_boson_operator`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ee6f8768ef75"
   },
   "outputs": [],
   "source": [
    "from openfermion.transforms import get_boson_operator, get_quad_operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "609b6483153d"
   },
   "outputs": [],
   "source": [
    "H = QuadOperator('p0 q0', 1/2.)\n",
    "H += QuadOperator('q0 p0', 1/2.)\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "30899a713341"
   },
   "outputs": [],
   "source": [
    "get_boson_operator(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "95ed065cfb41"
   },
   "source": [
    "Note that, since these conversions are dependent on the value of $\\hbar$ chosen, both accept a `hbar` keyword argument. As before, if not specified, the default value of $\\hbar$ is `hbar=1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "786ed58cea67"
   },
   "outputs": [],
   "source": [
    "H = BosonOperator('0 0^')\n",
    "normal_ordered(get_quad_operator(H, hbar=0.5), hbar=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "10a2e94ff205"
   },
   "source": [
    "## Weyl quantization and symmetric ordering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "02defd458960"
   },
   "source": [
    "We also provide support for the Weyl quantization - this maps a polynomial function of the form\n",
    "\n",
    "$$f(q_0,\\dots,q_{N-1},p_0\\dots,p_{N-1})=q_0^{m_0}\\cdots q_{N-1}^{m_{N-1}} p_0^{m_0}\\cdots p_{N-1}^{m_{N-1}}$$\n",
    "\n",
    "on the phase space to the corresponding combination of quadrature operators $\\hat{q}$ and $\\hat{p}$. To do so, we make use of the McCoy formula,\n",
    "\n",
    "$$q^m p^n \\rightarrow \\frac{1}{2^n} \\sum_{r=0}^{n} \\binom{n}{r} q^r p^m q^{n-r}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "39041da8098f"
   },
   "outputs": [],
   "source": [
    "from openfermion.transforms import weyl_polynomial_quantization, symmetric_ordering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3c24682b7ad5"
   },
   "source": [
    "For `weyl_polynomial_quantization`, the polynomial function in the phase space is provided in the form of a string,    where 'q' or 'p' is the phase space quadrature variable, the integer directly following is the mode it is with respect to, and '^2' is the polynomial power. If the power is not provided, it is assumed to be '^1'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "045dc06421f4"
   },
   "outputs": [],
   "source": [
    "weyl_polynomial_quantization('q0 p0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "066f8195101b"
   },
   "outputs": [],
   "source": [
    "weyl_polynomial_quantization('q0^2 p0^3 q1^3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "f19b276cc8d1"
   },
   "source": [
    "McCoy's formula is also used to provide a function that returns the symmetric ordering of a `BosonOperator` or `QuadOperator`, $S(\\hat{O})$. Note that $S(\\hat{O})\\neq \\hat{O}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "d7ae23880302"
   },
   "outputs": [],
   "source": [
    "symmetric_ordering(QuadOperator('q0 p0'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ea7958e22c94"
   },
   "source": [
    "Consider the symmetric ordering of the square of the bosonic number operator, $\\hat{n} = \\hat{b}^\\dagger \\hat{b}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "36130eb01dad"
   },
   "outputs": [],
   "source": [
    "from openfermion.hamiltonians import number_operator\n",
    "n2 = number_operator(1, parity=1) * number_operator(1, parity=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "d7c86afc29c7"
   },
   "outputs": [],
   "source": [
    "n2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "d3a2db524be9"
   },
   "outputs": [],
   "source": [
    "Sn2 = symmetric_ordering(n2)\n",
    "Sn2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2bbd7cdf62a3"
   },
   "source": [
    "We can use `normal_ordered_boson` to simplify this result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "3858191c9ea9"
   },
   "outputs": [],
   "source": [
    "Sn2 = normal_ordered(Sn2)\n",
    "Sn2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "50fb5cacb861"
   },
   "source": [
    "Therefore $S(\\hat{n}) = \\hat{b}^\\dagger \\hat{b}^\\dagger \\hat{b}\\hat{b} + 2\\hat{b}^\\dagger \\hat{b} + 0.5$. This is equivalent to $\\hat{n}^2+\\hat{n}+0.5$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7879325f20a2"
   },
   "outputs": [],
   "source": [
    "Sn2 == normal_ordered(n2 + number_operator(1, parity=1) + 0.5*BosonOperator.identity())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "172bf67ca8f6"
   },
   "source": [
    "## Bose-Hubbard Hamiltonian\n",
    "\n",
    "In addition to the bosonic operators discussed above, we also provide Bosonic Hamiltonians that describe specific models. The Bose-Hubbard Hamiltonian over a discrete lattice or grid described by nodes $V=\\{0,1,\\dots,N-1\\}$ is described by:\n",
    "\n",
    "$$H = - t \\sum_{\\langle i, j \\rangle} b_i^\\dagger b_{j + 1}\n",
    "         + \\frac{U}{2} \\sum_{k=1}^{N-1} b_k^\\dagger b_k (b_k^\\dagger b_k - 1)\n",
    "         - \\mu \\sum_{k=1}^N b_k^\\dagger b_k\n",
    "         + V \\sum_{\\langle i, j \\rangle} b_i^\\dagger b_i b_j^\\dagger b_j.$$\n",
    "         \n",
    "where\n",
    "\n",
    "- The indices $\\langle i, j \\rangle$ run over pairs\n",
    "  $i$ and $j$ of adjacenct nodes (nodes that are connected) in the grid\n",
    "- $t$ is the tunneling amplitude\n",
    "- $U$ is the on-site interaction potential\n",
    "- $\\mu$ is the chemical potential\n",
    "- $V$ is the dipole or nearest-neighbour interaction potential\n",
    "\n",
    "The Bose-Hubbard Hamiltonian function provided in OpenFermion models a Bose-Hubbard model on a two-dimensional grid, with dimensions given by `[x_dimension, y_dimension]`. It has the form\n",
    "\n",
    "```python\n",
    "bose_hubbard(x_dimension, y_dimension, tunneling, interaction,\n",
    "                 chemical_potential=0., dipole=0., periodic=True)\n",
    "```\n",
    "\n",
    "where\n",
    "\n",
    "- `x_dimension` (int): The width of the grid.\n",
    "- `y_dimension` (int): The height of the grid.\n",
    "- `tunneling` (float): The tunneling amplitude $t$.\n",
    "- `interaction` (float): The attractive local interaction $U$.\n",
    "- `chemical_potential` (float, optional): The chemical potential $\\mu$ at each site. Default value is 0.\n",
    "- `periodic` (bool, optional): If True, add periodic boundary conditions. Default is True.\n",
    "- `dipole` (float): The attractive dipole interaction strength $V$.\n",
    "\n",
    "Below is an example of a Bose-Hubbard Hamiltonian constructed in OpenFermion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2699f571f408"
   },
   "outputs": [],
   "source": [
    "from openfermion.hamiltonians import bose_hubbard, fermi_hubbard\n",
    "bose_hubbard(2, 2, 1, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "d58d4ec176e2"
   },
   "source": [
    "## Sparse bosonic operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "04673b1d4a9c"
   },
   "source": [
    "Like the fermionic operators, OpenFermion contains the capability to represent bosonic operators as a sparse matrix (`sparse.csc_matrix`). However, as the fermionic operators can be represented as finite matrices, this is not the case of bosonic systems, as they inhabit a infinite-dimensional Fock space. Instead, a integer truncation value $N$ need to be provided - the returned sparse operator will be of size $N^{M}\\times N^{M}$, where $M$ is the number of modes in the system, and acts on the truncated Fock basis $\\{\\left|{0}\\right\\rangle, \\left|{1}\\right\\rangle, \\dots, \\left|{N-1}\\right\\rangle\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "c8f604fa519b"
   },
   "outputs": [],
   "source": [
    "from openfermion.linalg import boson_operator_sparse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4a4895d9de71"
   },
   "source": [
    "The function `boson_operator_sparse` acts on both `BosonOperator`s and `QuadOperator`s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "a081c187c79a"
   },
   "outputs": [],
   "source": [
    "H = boson_operator_sparse(BosonOperator('0^ 0'), 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ab0fd120866b"
   },
   "outputs": [],
   "source": [
    "H.toarray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "54151b8cd511"
   },
   "outputs": [],
   "source": [
    "H = boson_operator_sparse(QuadOperator('q0'), 5, hbar=1)\n",
    "H.toarray()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "bosonic_operators.ipynb",
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
